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Abstract 

This paper is the first in a series that aims to understand the long-term 
evolution of neutron star magnetic fields. We model the stellar matter as 
an electrically neutral and lightly-ionized plasma composed of three moving 
particle species: neutrons, protons, and electrons; these species can be con- 
verted into each other by weak interactions (beta decays), suffer binary col- 
lisions, and be affected by each other's macroscopic electromagnetic fields. 
Since the evolution of the magnetic field occurs over thousands of years or 
more, compared to dynamical timescales (sound and Alfven) of milliseconds 
to seconds, we use a slow-motion approximation in which we neglect the in- 
ertial terms in the equations of motion for the particles. This approximation 
leads to three nonlinear partial-differential equations describing the evolu- 
tion of the magnetic field, as well as the movement of two fluids: the charged 
particles (protons and electrons) and the neutrons. These equations are first 
rather than second order in time (involving the velocities of the three species 
but not their accelerations). In this paper, we restrict ourselves to a one- 
dimensional geometry in which the magnetic field points in one Cartesian 
direction, but varies only along an orthogonal direction. We study the evolu- 
tion of the system in three different ways: (i) estimating timescales directly 
from the equations, guided by physical intuition; (ii) a normal-mode analysis 
in the limit of a nearly uniform system; and (iii) a finite-difference numeri- 
cal integration of the full set of nonlinear partial-differential equations. We 



1 



find good agreement between our analytical normal-mode solutions and the 
numerical simulations. We show that the magnetic field and the particles 
evolve through successive quasi-equilibrium states, on timescales that can be 
understood by physical arguments. Depending on parameter values, the mag- 
netic field can evolve by ohmic diffusion or by ambipolar diffusion, the latter 
being limited either by interparticle collisions or by relaxation to chemical 
quasi-equilibrium through beta decays. The numerical simulations are fur- 
ther validated by verifying that they satisfy the known conservation laws in 
highly nonlinear situations. 

1 Introduction 

Observations of the surface magnetic fields of neutron stars have shown a corre- 
lation between the age of the star and the magnetic-field strength. For instance, 
young radio pulsars and high-mass X-ray binaries have surface magnetic fields of 
the order of lO'^ G, while millisecond pulsars and low-mass X-ray binaries, which 
are older objects, have magnetic fields ~ 10^ G. These observations suggest that 
the magnetic field is decaying, although this might be a by-product of accretion 
(ID; 121; lEl; t4\). On the other hand, it is thought that the spontaneous decay of 
the ultrastrong 10^"*"^^ G magnetic field in magnetars is the main source of their X- 
ray luminosity since these objects appear to radiate substantially more power than 
that available from their rotational energy loss (|5 |; [6|; [7|). Since these objects 
appear to be isolated, the field decay must be attributed to processes intrinsic to the 
neutron stars. 

Ref . m ; (hereafter GR-92) identified three long-term mechanisms that can pro- 
mote the spontaneous decay of the magnetic field in the interior of neutron stars: 

1. Ambipolar diffusion, i.e. the drift of the magnetic field and the charged par- 
ticles relative to the neutrons. Ref. |6|; (hereafter TD-96), using magnetar 
parameters, estimated the timescale of the magnetic field decay by means of 
ambipolar diffusion, finding that was consistent with the timescale of 10^~^ 
years observed for these objects. 

2. Hall drift, a non-dissipative advection of the magnetic field by the associated 
electrical current (lEl; m; IE]; 113, O; 111; IBl; III, lEl; ITSl). 

3. Ohmic diffusion, a dissipative process caused by electrical resistivity (see 
also Baym et al. Il9l ). 

All of these processes occur over timescales of thousands of years or more, com- 
pared to typical dynamical timescales (sound and Alfven) of miliseconds to sec- 
onds. The work of GR-92 was analytical, and therefore useful to identify general 
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processes and relevant timescales, but it did not address the action of the identi- 
fied processes in their full nonlinear development and their interactions with each 
other. The full evolution of the magnetic field can only be addressed by numerical 
simulations. 

Regarding the geometrical configuration of the magnetic field, three-dimensional 
magnetohydrodynamics (MHD) simulations showed that in a stable stratified, non- 
rotating star, the initial field evolves on a short, Alfven-like timescale to a large- 
scale equilibrium configuration that can be axisymmetric ( 112^ . ETI : |[22l ) or non- 
axisymmetric (|23l), depending on the radial dependence of the initial magnetic- 
field strength. These simulations, being focused on the evolution of the magnetic 
field over short dynamical timescales, were carried out within the framework of 
the MHD theory, which considers the stellar matter to be a single fluid. The evo- 
lution of these configurations has not yet been studied over longer timescales, on 
which relative motions of different particle species are present and the processes 
studied by GR-92 become important. The description of these long-term processes 
requires a multi-fluid theory. 

In this paper we extend the model of GR-92 by allowing both neutron and 
charged-particle movements, rather than considering the neutrons as a fixed back- 
ground (see Sect.©. For simplicity, we restrict ourselves to a plane-parallel system 
where the magnetic field points in one Cartesian direction but varies only along an 
orthogonal direction. Although this one-dimensional model is unrealistic, it allows 
us to study the basic evolutionary processes in an analytical and numerical fash- 
ion with no serious complications before considering a realistic stellar geometry in 
more dimensions. GR-92 considered all the species in the star as non-interacting 
fermions; in contrast, we use equations of state including interactions among pro- 
tons and neutrons. Since in neutron stars the magnetic pressure is much less than 
the degeneracy pressure of the particles, we treat the particle densities in terms 
of small perturbations around a non-magnetized background in full equilibrium. 
Thus, we obtain a system of equations that is nonlinear with respect to the mag- 
netic field but linear with respect to the particle densities (see Sect.O. 

We will show that the magnetic field and particle species evolve through suc- 
cessive quasi-equilibrium states. We describe the basic physical processes of this 
evolution and estimate analytically the characteristic timescales required to reach 
these quasi-equilibrium states. We find that these timescales depend on the val- 
ues of physical parameters in the system, such as the magnetic-field intensity, the 
collision rate between the charged particles and the neutrons, the weak interaction 
rate, and the ohmic resistivity (Sect. 14.11 ). To obtain an analytic solution of the 
system of partial differential equations, we linearize the equations with respect to 
the magnetic field and find normal-mode solutions (Sect. 14.21 ) and obtain analytical 
approximations for the decay times in Appendix 16.11 We develop numerical sim- 
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ulations (whose basic algorithm is described in Appendix 16.21 ) in the hnear regime 
and compare with the normal-mode results (Sect. 14.31 ). In Sect. |4.4[ we carry out 
numerical simulations in the nonlinear regime and verify the conservation laws of 
our model. Finally in Sect. [5j we summarize the main results of this paper and 
provide our conclusions. 

2 General Physical Model 

We model the stellar interior as an electrically-neutral and slightly-ionized plasma 
composed of three moving-particle species: neutrons («), protons (p), and electrons 
(e). This is an extension of the model of GR-92, since these authors considered the 
neutrons to form a motionless background. Furthermore, we account for strong 
interactions between neutrons and protons by writing each of their chemical poten- 
tials as a function of both of their number densities: yu„,p(r, t) = fj.n,p[nnir, t), np{f, f)]. 
We consider the electrons to be an ideal, relativistic Fermi gas, and write their 
chemical potential as a function only of their own density, yUe(r, t) - fieirieif, t)) = 
yJimgC^)^ + ffic^{?>n^ne(f, 0]^^^, with nig the electron rest mass. The plasma species 
are described as three fluids coupled by collisions and electromagnetic forces, sat- 
isfying the equations of motion (GR-92; Reisenegger et al. |16]): 



where v,- is the mean velocity; djdt = d/dt + (v, • V), yU,/c^ is the effective mass 
of each species (which includes corrections due to interactions and relativistic ef- 
fects; [i24J ); rii^fii is the degeneracy-pressure gradient of the species /; ifr is the 
gravitational potential; E and B are the electric and magnetic fields; and the last 
term represents the drag forces due to elastic binary collisions, which damp the 
relative motions of the different species. The collisional coupling strengths are 
parametrized by the symmetric matrix y,-y that depends generally on the local den- 
sity and temperature. We ignore the effects of superfluidity or superconductivity, 
as well as additional particles that might be present in neutron stars. 

In the early stages of the neutron-star formation, particles are locked by colli- 
sions, which behave as a single fluid. Alfven waves propagate across the star, al- 
lowing it to reach a magnetohydrostatic quasi-equilibrium state in which all forces 
acting on a fluid element are closely in balance. Since in this paper we are inter- 
ested in simulating the evolution over far longer timescales characteristic of the 



rii — — = -riiVni -rii — Wi// 
dt 




(1) 
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processes of GR-92, we use a slow-motion approximation in which we neglect the 
acceleration terms on the left-hand side (LHS) of Eq. ([T|) (i.e., (?i,;u,/c^) dvj/dt = 0). 
We assume implicitly that all the forces acting on a given particle species balance 
each other at all times. To some extent, this will happen, in the sense that the rela- 
tive velocities will take such values that the drag forces balance all the other forces, 
and these drag forces will decay as the other forces come into balance. However, 
if Eq. ([T]) is summed over all three particle species, the drag forces cancel and 
no mechanism remains to ensure that the entire fluid reaches magnetohydrostatic 
quasi-equilibrium. To create such mechanism, we introduce an artificial friction- 
like force acting on the neutrons (the most abundant species) of the form -rinOVn, 
such that the neutron equation of motion becomes 



= -UnSUn - rin-^^i^ - ^ JnjrinnjiVn - Vj) - n„aVn. (2) 

As a result, we obtain the equation of magnetohydrostatic quasi-equilibrium for a 
fluid element. 



V A*! ,4 JxB 

niVni- 2_j nj — yif/+ n^avn, (3) 

^p,n i=e,p,n 



where J = nce{v*p - v^) - {c/4n){^ x B) is the electrical current density. Thus, the 
system approaches magnetohydrostatic quasi-equilibrium on a timescale controlled 
by the parameter a (see Sect. 14.1.11 ). This pai^ameter is chosen in such a way 
that this timescale is sufficiently long for the numerical code to be able to deal 
with it (and therefore much longer than the true dynamical timescales [sound and 
Alfven]), but shorter than the long timescales of interest in this paper. 

If we add Eq. ([T|l for electrons and protons, assuming charge neutrality, rig = 
Up = fic, we obtain the diffusion equation for the combined fluid of charged parti- 
cles, 

ncrinJcnVA ^ -ric" t^c - Tlc^Vlf/ + , (4) 

c-^ c 

where //^ = f^e+Mp, Jen = 7en + 7pn, ^nd we define the ambipolar diffusion velocity 
by: 

^ _ ypn(Vp - V„) -I- y,niVe - V„) 

Va = • (5) 

yen 

If the charged particles are not in diffusive quasi-equilibrium, that is the forces on 
the right-hand side (RHS) of Eq. ^ are not in balance, the latter will drive them to 



5 



move relative to the neutrons at the ambipolar diffusion velocity va- On the other 
hand, if the system is in magnetohydrostatic quasi-equilibrium (see Eq. Q and 
the charged particles have reached their diffusive quasi-equilibrium (va = 0), this 
implies that neutrons are also in a diffusive quasi-equilibrium state given by 

6 = -Un^Hn - tln^Vl//. (6) 

To derive the equation that governs the evolution of the magnetic field, we com- 
bine Eq. ([T]), for electrons and protons without the inertial terms, and the induction 
equation V x £" = -{\/c){dB/dt), to obtain 



— = Vx[(i/„ + i?A+Vff)xBj-Vx 

C^(7en-7pn\ ^ 
-—V X VyU, 

2e \ yen I 

I ^(yenf^p-7pnHe\ ^, 

V X V(A, 

ec \ 7cn I 



c^V X 
Ancr 



(V) 



where we define the Hall drift velocity, which is proportional to the electrical cur- 
rent density, by 



Vh 



7en 7pn 



C{7en - 7pn) 



(8) 



7cn nce7cn 
The first term on the RHS of Eq. Q shows that the magnetic field is transported 
by the sum of the neutron velocity Vn, the ambipolar diffusion velocity va, and the 
Hall drift velocity v//. The second term represents the ohmic diffusion, where the 
electrical conductivity is 



cr 



7pn 7en 



-1 



(9) 



Finally, the last two terms on the RHS of Eq. ([7]l represent battery effects. 

To complete a set of equations to describe the system evolution, we need equa- 
tions for the evolution of the particle densities. Weak interactions between the 
particles tend to reduce chemical-potential imbalances between the charged parti- 
cles and neutrons. We define the differences between the rates, per unit volume, of 
the reactions p + e ^ n + Ve and n ^ p + e + 'v^ as AT = T{p + e ^ n + Vg) - T{n — > 
p+e + v^) = /lA/i, where A// - fic-^^n is the chemical imbalance, and the parameter 
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A depends generally on density and temperature (GR-92) 0. The chemical quasi- 
equilibrium between the charged particles and neutrons is reached when A// - 0. 
The continuity equations for the particle densities are then given by 

^ + ^ ■{niv-'i) = ±A{^^l). (10) 

ot 

In Eq. ([Tol l, the + sign corresponds to the neutrons and the - sign to the electrons 
and protons. Defining the baryon number density as = n„ + ric and adding 
Eq. ([Tol l for charged particles and neutrons, we obtain the conservation law for ng, 

— + V • (rinv'n + ripVp) = 0. (11) 

The set of differential Eqs. (|7]l, dTOl ). and ([TTI ) are first, rather than second order 
in time, that is they involve the velocities of the three fluid species, but not their 
accelerations. 



3 One-Dimensional Model 
3.1 Basic Equations 

We consider a one-dimensional geometry in which the magnetic field points in one 
Cartesian direction z, but varies only along an orthogonal direction x as B{r, t) = 

By{x, t)z. We also assume that all of the other physical variables vary only along 

-* 

X, and therefore that the gradient operator is V = x{d/dx). From Ampere's law, 
Jjc = {c/4n){V X B)x ^ ^ riceivpx - Vex), thus, Vgx = Vpx = v„ = v„x + vax- 
Using Eqs. ©, Q, (W^, and ([Ull, and defining r = c^/iAncr), we obtain the 
following non-linear set of equations for the evolution of the magnetic field and 
particle densities: 

dB,_ divM ^d_(dB,\ ^^2) 



dt dx dx\ dx 

-r- = -T- (rinVnx + flcVcx) , (13) 

Ot ox 
drir d 

-^ = -—(n,v,x)-A(Aix), (14) 
ot ox 



where 



'If A/i > 5kT, where k is the Boltzmann constant and T is the temperature, A must also be allowed 
to depend on Ayu (Reisenegger 1251 : Fernandez & Reisenegger 1261 ). 
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an„ 



dx dx 



+ n„ 



dx dx 
dx\%n 



and 



(15) 



1 



VAX = — 



nnncjc. 



dx dx j dx Stt 



(16) 



In this one-dimensional geometry V xivn x Bj - 0, and, because we are consider- 
ing only variations along x, the gradients that appear in the battery terms of Eq. (|7]l 
are parallel. Thus, there is no contribution from the Hall drift and the battery terms 
to the evolution of the magnetic field. 



3.1.1 Conserved quantities and boundary conditions 

We note that Eqs. ([T2l ) and ( fT3] ) can be written as flux-conserving equations of the 
form 



dA dS, 



dt dx 

where A represents either B^otub, and 5 a is defined in each case by 



(17) 



dB 

^ B- - VcxBz - ■ 

and 



Sb. = VcxB,-—^, (18) 
47ro" dx 



Sns ^nnVnx + ricVcx. (19) 

On the other hand, if we define Ca = ^^"""^ A{x, t)dx, we obtain the conserva- 
tion law 

d 

— Ca - ^aIa-o - 5'aU=jc„„, (20) 

We note that Cb. = ^sit) is the magnetic flux, while C„g = A^b(0 is the baryon 
number. We see from Eq. (|20l ) that Ca is a conserved quantity if S a\x=o-S a\x=x„„, - 
0, which depends on the boundary conditions imposed on the variables contained 
in 5 A- To conserve both the magnetic flux and the baryon number during evolution, 
we impose the boundary conditions 
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Vcx(x = 0, = Vcx(x,ruix, t) = 0, 



(21) 



Vnx(x = 0, ?) = VnxiXmax, t) = 0, (22) 

—^{x = 0, = -^{x = t) = 0. (23) 
ox ox 

Of course, this restriction does not provide a realistic description of a neutron star, 
in which magnetic flux can be lost through the boundaries; however, it allows us 
to control the precision of the calculation in this one-dimensional case. Although 
this geometry is quite unrealistic, it enables us to study the timescales involved in 
an analytical and numerical fashion. In the future, when we study this system in 
higher numbers of dimensions, we will use more appropiate boundary conditions 
for the magnetic field. 



3.2 Linearization of the equations with respect to the densities 



Since in neutron star conditions, the ratio between the magnetic pressure B-^/Sn 
and the degeneracy pressure of the charged particles is very small, we consider 
the magnetic field to be a small perturbation on a non-magnetized background in 
full hydrostatic, chemical, and diffusive quasi-equilibrium. If we label the phys- 
ical variables characterizing the background by the sub-index 0, the hydrostatic, 
diffusive, and chemical equilibria are given by 

^.^f^ = 0, (24) 
ox ox 

ox ox 

A/lo = /lOc - y"On = 0, (26) 

The magnetic field cannot force significant displacements of the particles; we 
therefore write ?i,(x, i) = rn)i{x) + 6ni{x, t), where \6ni{x, t)\ «; noiix). The chemical 
potentials are fii(x, t) = ^Qi{x) + 6ni{x, t), where 5ni{x, t) - kiBSrisix, O+kicdndx, f), 
kiB - (dfii/driB),^^, and kic = (5yU,/5?ic)„og, with keB = and the remaining coeffi- 
cients are calculated from the equation of state. We also use the Cowling approx- 
imation, neglecting the perturbations of the gravitational potential with respect to 
the background value, i.e., ifr = iJ/o.By neglecting terms of second order in Jn,- and 
using that Vcx - ^nx + ^Ax, ^ns = (5n„ -l- dric, and the definitions % = kcc - Kc, 
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kcc = kpc + kec, koB = kpB - k,,^, and ro = c /(Attcto), we derive the following set of 
equations, which are nonlinear with respect to the magnetic field By{x, t), but linear 
with respect to the density perturbations <5«b(x, t) and dndx, t), namely, 

+ T-\ro—\, (27) 



dt dx dx\ dx 

dduB d 

— r— = -T- {riQnVnx + "OcVcx) , (28) 

ot ox 
ddric d 

-H— = (nocVcx) - ^ ihocdric + koBSriB) , (29) 
ot ox 



where 



anon 



nOnl^On-^ 

ox 



d ( KJub + kncSric 



d (kccSric + kpB6nB\ d (Bj 



(30) 



and 



VAX 



nonnocTcn 



"Oc/^OnT- 

ox 



d I kccSric + kpBduB 



.mi 

dx[S7r) 



(31) 



3.3 Dimensionless equations 

We proceed in writing the set of Eqs. ([27]), (l28]l, ([301), and ([SB in terms of 
dimensionless variables. We start by defining dimensionless variables as a = a/a.v> 
with the sub-index s that represents a characteristic value of the corresponding 
dimensional variable. In the following, we explain the meaning of each of the 
characteristic values. We normalize the space variable x with respect to the length 
of the system, that is the length of our computational domain d, thus x^ = d. We 
note, however, that there is another characteristic length L, which is the length 
over which the functions vary spatially; we therefore use this scale when deal- 
ing with order-of-magnitude estimates of spatial derivatives. The time variable t 
is normalized with respect to tg = ad^ /{nonk„B). The meaning of this quantity 
will become clearer in the following section; it is related to the shortest relevant 
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timescale of the system. We assume that Bg - B"'"^ is the maximum value of 
the initial magnetic field, ris = is the maximum neutron background den- 
sity, rics = «Q is the maximum of the charged background density, = k'"^'' and 
= yu™^ are properties of the background, = r[J'"^ is the maximum resistivity in 
the background. From the third term on the RHS of Eq. we estimate a charac- 
teristic order-of-magnitude value for the neutron velocity induced by the magnetic 
stress to be - B^JiSnn^Xsa). If we compare the first and third parentheses on 
the RHS of Eq. (I30l ). and use that close to magnetohydrostatic quasi-equilibrium 
dng/noB ~ dric/riQc, "on/'^OB ~ !> ^rid nQ„k„B » n^^kn^, then we can estimate a char- 
acteristic neutron-density perturbation to be 5ns - B^J{^jmsks). Finally, defining 
ncs = ncs/ris and k = k/k^ where k is any of the corresponding parameters defined 
in the last section, we obtain the following dimensionless set of equations. 



where 



^^-n T^(Z££^ + CO— |?o— I 
dl dx dx\ dx ' 

dSns 9 - - - X 



dSRc 



dl 51 



= - — {nocVcx) - ^csO (kocSEc + koBSns) , 



(32) 
(33) 
(34) 



1 



_ _ d ( knsSnB + kncSric 
dx { 



d 

dx 



2 

' kccSEc + kpBSriB \ d(B, ) 
dx 



(35) 



VAX 



d {kccSEc + kpBSns^ 



dx 



(36) 



In the above set of equations, the parameter Y = B^/(87r«ocj«o«4^;iBj) controls 
the coupling strength between the magnetic field and the particle species, ca = 
arg /(risks) controls the importance of the ohmic diffusion, = aAgX^l^cs controls 
the weak interaction rate, and finally e = alinsjs) controls the importance of the 
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ambipolar diffusion process. In diffusive quasi-equilibrium, Eq. (I30l ) witli Vnx ~ 0, 
and Eq. ( [3T| ) witli vax ~ 0, imply that 

— = T, (37) 

where we used that the fraction \nonknB\/\nocikcc + kps)] ~ 1. For instance, in 
the non-interacting particle limit of GR-92 this fraction is approximately equal 
to 2 since the chemical potential of each particle species / is a function only of 
its number density t) = fii{ni{f, t)) = yj{miC^)'^ + ff-c\?>n'^ni{r, 1)^1^, which 
implies that knB = {diUn/dn„)^^^^, knc = -kns, kpB = 0, and kj,. = {dmldric)^^^,. On 
the other hand, for the background numerical values that we use in this paper this 
fraction is w 0.7. These background parameters are extracted from the equation 
of state provided by Akmal et al. (|24|), evaluated at a neutron number density 
riQn - 9.8 X 10^^ cm~^, which corresponds to a charged-particle number density 
of wqc - 3.9 X IQ^^cnT^ and to a total mass density (including neutrons, protons 
and electrons) of 1.5 x 10^"* g/cm^. For these values, k^B ~ 4.2 x lO'^^erg cnv', 
knc w -1.2 X IQ-'^^erg cnr' , kpB « -8.2 x IQ-'^^erg cm^, kp^ « 2.8 x IQ-'^^erg cwr', 
and kec w 1.3 x 10""*^ erg cw?. In addition, the linearization with respect to the 
variable dric requires that \6nc\ «c noo i-e-> T «; 1 and therefore 6.3 x lO^^G. 



4 Results 

4.1 Characteristic evolutionary timescales 

We consider a dynamical system of three independent variables controlled by three 
differential equations that are first order in time. Thus, in the linear limit, the 
system will have three exponentially decaying "modes" of different timescales. 
Although, this is not strictly true in the general nonlinear case, one can always 
identify three characteristic timescales on which the system approaches successive 
quasi-equilibrium states. Real neutron stars approach a magnetohydrostatic quasi- 
equilibrium on a timescale not much longer than the Alfven time, which is of the 
order of seconds. Here, as described in Sect. |2l this early evolution is mimicked 
by a artificial friction force term proportional to a parameter a; this parameter is 
chosen so that this timescale is sufficiently long enough for the numerical code to 
be able to deal with it, but shorter than the timescales of modeled processes, such 
as ambipolar diffusion and weak interactions. In what follows, we make analytic 
estimates of the characteristic evolutionary timescales and evaluate them for typical 
magnetar core parameters (see e.g.. Arras et al. [7|). In doing this, we can make 
sense of our numerical results and estimate the order of magnitude of the timescales 
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involved in a real system, under the limitations of our one-dimensional model. 



4.1.1 Timescale to achieve magnetohydrostatic quasi-equilibrium 

We assume that our non-magnetized background star is in hydrostatic quasi-equilibrium 
[see Eqs. (I24l ) and (|25l)1. in other words, the net force on a fluid element containing 
all species is zero. When a magnetic field is present in the system, the magnetic 
force applies pressure to the charged particles (electrons and protons), inducing 
density perturbations that create an imbalance between the diff'erent forces acting 
on a fluid element [see Eq. (l30l)1. The particles must move until a magnetohy- 
drodystatic (MHS) quasi-equilibrium is reached. The neutron velocity necessary 
to achieve this balance is expressed in Eq. (I30l) . During the first stages of evolution, 
the collisional coupling between the charged particles and neutrons compels them 
move with about the same velocity, Vcx ~ Vnx- Neglecting weak interactions, from 
Eqs. |33]and|34]we obtain the consistency condition 



SuB/noB » Sric/noc 1- (38) 

The induced charged-particle and neutron-pressure gradients tend to choke the 
magnetic force. The magnetohydrostatic quasi-equilibrium state is reached when 
there is a close balance between these opposing forces [see RHS of Eq. (I30l)1. that 
is. 



non I knB + —Kc W«B ~ , (39) 

\ riQB I Sn 

where we neglected the second term in the RHS of Eq. ( |30l ) since tiqc non 
and used the Eq. ([38] ). The velocity induced by the initially-unbalanced magnetic- 
pressure gradient is 

Vnx ~ Vcx ~ Z '~~r- (40) 

On the other hand, from Eq. (|28]) . we obtain the time required to create (or destroy) 
a perturbation Shb as 

LSnn 

tMHS -■ (41) 

nOBVnx 

Using Eqs. (l39l ). (l40l) . and ((4TI) . and uob ~ «o« we estimate the timescale to reach 
the magnetohydrostatic quasi-equilibrum as 

tMHS 7 — « — 1— - b ' (42) 
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thus, Imhs ~ (L/d)^ < 1. In Eq. (I42l ) we used that |(«oc^nc)/(«On^nB)l ^ 1- In fact, 
in the non-interacting particle limit of GR-92 we have k„c = -Kb, thus, this frac- 
tion is of the order of |?ioc /"OmI ^ 1- For the numerical background values that we 
use in this paper, this fraction is ~ 0.12. We note that the scale Imhs is the short- 
est relevant timescale in the system, and is controlled by the artificial a parameter 
that was introduced in our slow -motion approximation (Sect.O. Real neutron stars 
evolve to the magnetohydrostatic quasi-equilibrium on a short timescale not much 
longer than the Alfven time, which for typical magnetar core parameters scales as 

tAlfven = 5.1 X 10'^ B\l, S, (43) 

where = R/(10^ cm) denotes the radius of the star in units of 10^ cm, and 
Bi5 = B,/(10^^ G) the magnetic field in units of lO'^ G (this timescale as well 
as the following ones are evaluated at the typical mass density 1.5 x 10^^ g/cm^). 
This timescale is far shorter than the timescales of the processes that promote the 
long-term evolution of the magnetic field (see Sects. I4.1.2H4.13] ). which are of the 
order of years or much longer. A numerical code simulating the evolution on the 
Alfven timescale would require a time step many orders of magnitude shorter than 
that required to simulate the long-term evolution in a computational time that is not 
prohibitively long. We overcome this difficulty by replacing the short-term dynam- 
ics by the artificial friction term proportional to the parameter a. This parameter is 
chosen so that the timescale in Eq. (I42l ) is long enough for the numerical code to be 
able to deal with it (and therefore much longer than the Alfven time), but shorter 
than the timescales of the long-term processes that we discuss in the following 
sections. 



4.1.2 Timescale for charged particles to reach diffusive quasi-equilibrium 
through ambipolar diffusion neglecting weak interactions 

We now assume that the magnetohydrostatic quasi-equilibrium, discussed in the 
last section, has been reached. However, the charged particles continue to move rel- 
ative to the neutrons due to ambipolar diffusion (subject to collisional drag among 
different species), and they reach diffusive quasi-equilibrium when there is a close 
balance between the magnetic force and the charged-particle pressure gradients 
(see RHS of Eq. (|3TI ). By an argument analogous to that in the previous section and 
using the temperature dependence of the collisional frequencies from Ref. 1271 . we 
obtain the timescale tdrag for charged particles to reach diffusive quasi-equilibrium 
as 
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where L5 = L/(10^ cm) and Tg = T/{10^K). We can write this time in units of our 
normahzation time as 

- /L\2 no„ycn ~tMHS /^\^ 1 

''"■"^ ~ I J j ^ r ~ I J j ?' ^^^^ 

We require Idrag ^ ^m//s> therefore ? «; 1. The ? parameter in Eq. (1451 ). which 
is inversely proportional to the collisional parameter y„i> controls the timescale on 
which the charged particles reach the diffusive quasi-equilibrium. 



4.1.3 Timescale to achieve chemical quasi-equilibrium through weak inter- 
actions, neglecting ambipolar diffusion 

If there is a perturbation of the chemical quasi-equilibrium (kocdric 4^ -koB6nB; 
see RHS of Eq. (|29l ). the characteristic timescale on which the chemical quasi- 
equilibrium is restored through weak interactions (charged particles decaying into 
neutrons and viceversa) can be estimated from Eq. ( [291 ). Neglecting the first term 
on the RHS, which takes into account the ambipolar diffusion, comparing the terms 
with dric, and using the temperature dependence of the A parameter from Ref. ESl . 
which assumes modified Urea reactions (e.g., Ref. [29]), we obtain 



tweak 



1 



4.3 X 10^ T-^ yr. 



A{koc + koB) 

If we write Eq. (l46l ) in units of our normalization time, we derive 



tweak 



nOnknB 



aAcP- 



1 



(46) 



(47) 



aAd^koc + koB) 

Using Eq. (|47] ). we obtain Iweakl'tMHS ~ id/L)-^{l/9); we therefore require 9^1, 
so that tweak ^ tMHS ■ The 9 parameter in Eq. (|47] ) is directly proportional to 
the weak interaction rate parameter A, which controls the timescale on which the 
chemical quasi-equilibrium is achieved. 

Both weak interactions and ambipolar diffusion processes contribute in general 
to the decay of the charged-particle density perturbations. The more rapidly-acting 
processes determine the evolutionary timescale of the charged particles. We pro- 
pose an approximate general interpolation formula: 



tSllc 



1 1 


-1 




z 1- z 






. tdrag tiveak . 







(48) 
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4.1.4 Ohmic diffusion timescale 



The timescale on which the magnetic field decays by ohmic diffusion can be esti- 
mated from Eq. dTTl ). neglecting the first term on the RHS and using the temperature 
dependence of the electrical conductivity from Ref. QUI , as 

tohmic ~ — = ^^^^ ~ 1.4 X 10" Lj r,-2 yr. (49) 
On the other hand, using Eq. (B^ . we obtain 

/L\2 nonknB ~tMHS (W ^ r<c^^ 
tohmic ~ "7 = "7 = (50) 

\ai aro to \a/ a> 

The a» parameter is directly proportional to the parameter ro and controls the ohmic 
diffusion timescale. We require 73 <^ 1, so that lohmk ^ (mhs- 



4.1.5 Timescale for the evolution of the magnetic field by ambipolar diffusion 
with weak interactions 

The two terms on the RHS of Eq. ( [32l) set the two basic timescales on which the 
magnetic field decays. The first of these terms couples the magnetic-field evolution 
with the particle dynamics, while the second one gives the magnetic-field evolution 
due to ohmic diffusion. For the remainder of this section, we estimate the magnetic- 
field evolution timescale by neglecting the second term (formally we set a» = or 
'tohmic o°)- Again, it is useful to distinguish between the two opposite extreme 
regimes discussed in Sects. |4~L2] and 14.1.31 

In the first regime, the ambipolar diffusion occurs more rapidly than the weak- 
interaction processes, that is tMHS t^irag t^eak- First, the system reaches the 
magnetohydrostatic quasi-equilibrium in the short timescale tMHS ■ In the second 
stage, all particles reach the diffusive quasi-equilibrium in the timescale t^rag and 
there is a close balance between the Lorentz force and the charged particle pressure 
gradient [see Eq. (|3TI)1. Since the magnetic field generates density perturbations in 
the charged particles but not in the neutrons, it prevents the system from reaching 
chemical quasi-equilibrium. Weak interactions tend to restore the local chemi- 
cal quasi-equilibrium on a characteristic timescale t^eak- This tendency causes a 
shght reduction of the pressure gradient in the charged particles with respect to 
the Lorentz force, compelling the charged particles and the magnetic flux to move 
together at a small ambipolar velocity vax ~ v^. This movement maintains the 
charged-particle pressure gradients, which tend to be erased by weak interactions 
that continue to operate. This interplay of weak interactions and ambipolar diffu- 
sion stops only once the pressure and magnetic field gradients disappear. 
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To estimate the timescale of this process, we first note that the ambipolar ve- 
locity needed to keep the charged-particle density perturbations stationary with 
respect to time {dSnddi - 0) can be estimated from Eqs. ( |29l ) and (l37l) as 

LBjAikoc + koB) 
VAX ^ ; • (51) 

oJmOcnOnKnB 

Since we neglect ohmic diffusion, the magnetic-field evolution is governed by the 
coupling with the charged-particle movement at the ambipolar velocity given by 
Eq. (ISTI ). From Eq. (|27] ). and using Eq. (|5TI ). we can estimate the timescale on 
which the magnetic field evolves in this case to be 

(1) ^ ^weak 



ambip Y BU{koc+koB) 

~ 1.7x lO'' T-^ yr. (52) 
Using Eqs. (|37] ) and (|52] ). we write this timescale in dimensionless form as 

\imhip ~ ^^^^ 

In this regime, the weak interactions operate very slowly, converting particles 
of one species into the other (beta decays) in a tendency to erase the chemical 
imbalance. The beta decays perturb the diffusive quasi-equilibrium, promoting 
magnetic-field evolution by means of ambipolar diffusion in a timescale that is 
limited by how rapidly the weak interactions operate (tweak)- We note that the 
timescale given by Eq. (l52l ) is the timescale to reach chemical quasi-equilibrium, 
tweak, amplified by the factor 1 /T (which is of the order of the ratio of the charged- 
fluid pressure to the magnetic pressure). The dependence of the timescale on this 
factor is expected since the Lorentz force drives the ambipolar diffusion, which 
maintains the chemical imbalance as long as there is a magnetic-field gradient. 

In the opposite regime, the weak interactions occur much faster than the am- 
bipolar diffusion process, i.e. tMHS tweak t^rag- As before, during the first 
stage of this regime, the system reaches the magnetohydrostatic quasi-equilibrium 
in the short timescale tMHS ■ At the end of this stage, the gradient of the magnetic 
pressure is balanced by the combined degeneracy pressure of all particles, with the 
neutrons providing the main contribution due to the much higher density. During 
the second stage, chemical quasi-equilibrium is established at every point in the 
system, coupling the neutron and charged-particle density perturbations. Since the 
charged-particle pressure is far smaller than that of the neutrons, it can not by itself 
balance the magnetic -pressure gradient, which, according to Eq. (|3TI ). causes an 
ambipolar drift of velocity 
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On the other hand, from Eq. (|27] ) we have tambip ~ L/vax, thus, we can estimate 
the magnetic-field evolutionary timescale in this case to be 



ambip 



L 



~ 1.8x10^ BllL]Tlyr. 



(55) 



If we use Eq. (1451 ). we obtain 



ambip 



( 




(56) 



In this case, the beta decays restore the chemical quasi-equilibrium quickly 
but the collisions between particles are frequent. This prevents the particles from 
moving easily relative to each other to achieve diffusive quasi-equilibrium. In this 
regime, the main mechanism promoting magnetic-field evolution is ambipolar dif- 
fusion. This process carries the magnetic field with a velocity that is limited by the 
collision rate between particles and depends on the magnetic-field strength. Equa- 
tion. (|55] ) gives the magnetic-field evolutionary timescale in this regime, which is 
that required for the particles to reach diffusive quasi-equilibrium, tdrag, but, as 
above, amplified by the factor 1/T. 

We conclude that, in general, where both weak interations and ambipolar dif- 
fusion contribute to the magnetic-field evolution, the slower of these processes de- 
termines the evolutionary timescale of the magnetic field. For a general estimate, 
we use the approximate interpolation formula that recovers the timescales of the 
regimes discussed above to the corresponding limits 



Equation (|57] ) corresponds to the estimate of Eq. (35) in GR-92 in the limit of non- 
interacting particles. For a case that includes ohmic diffusion, the shorter of the 
timescales given in, Eqs. ( |49l ) and (|57] ) provides the timescale for magnetic-field 
evolution, which we estimate using the interpolation formula 



tambip ^ ambip ^ ambip' 



(57) 



1 1 



1-1 



t ambip ^ohmic 

that can be written in dimensionless form as 




(58) 
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tB ~ 



1 1 

J/ 6 



n-1 



(59) 



The ordering between I^,,^ and 1b depends on the relative size of To and the other 
dimensionless parameters. We can place the relevant timescales in increasing order 
as 



t\ - tMHS, 



(60) 



t2 ^ min{f5„ . ffi) w I J- + ^ 



(61) 



(62) 



where 7i «: ?2 < ^3. It is easy to show that, for oj = 0, the ordering between the 
timescales is ?i = 1mhs> h = hn^' and ?3 = 1b- 

4.2 Normal Modes 

The set of Eqs. (l32l ). (l33l) . and (|34] | is nonlinear with respect to the magnetic field 
variable B^. To find a linear solution of this set, we assume that the properties of 
the background star are homogeneous, i.e., quantities with sub-index zero do not 
depend on position, and we linearize the magnetic field as 



B^{x, t) ^ B, + 6B^{x, t). 



(63) 



where Be is a constant and \6B^{x, t)\ is a small magnetic perturbation. If we 
choose Be as the magnetic-field normalization unit (i.e., B, - B^/Bc), we obtain 



B^CxJ) ^ \ + 6Bz(x,l). 



(64) 



We find normal-mode solutions for the dynamical variables in this linear sys- 
tem as 



SnB ^ a^g(??^)cos(/7rx)exp(-77,„?). 



(65) 



(5«c = «Si,('7m)cos(/;rJ)exp(-77„0, 



<5Sz = a^B i^m) cos(lnx) exp(-?/ J), 



(66) 
(67) 
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where m - 1 , 2, 3 is an integer labeling each mode, t„, - 1 /?7^ is the decay time 
of each mode, I is an integer measuring the mode wave-number, and agp is the 
amplitude of the corresponding variable 5p in each mode. We note that the func- 
tion cos(/7r^) satisfies the boundary conditions imposed in the previous section and 
that the characteristic length for spatial variations in this normal-mode solution 
is L = d/{ln). We note that the linear decay times fm correspond to the esti- 
mated characteristic timescales in the general nonlinear regime [see Eqs. (l60l ). 
(|6TI ). and (l62l)1. and we use the same notation in the remainder of this paper. The 
normal-mode solution translates into an eigenvalue problem that can be solved ei- 
ther numerically or using analytical approximations. In Figs. \VlQ\ we compare the 
decay times of the normal modes obtained from an analytical approximation to the 
eigenvalue problem (open triangles; see Appendix 16. II ) with those obtained from a 
numerical solution of the eigenvalue problem (black diamonds). We observe that 
in most cases there is good agreement between these two methods. Figures. [Ub) 
and [TJd) show a discrepancy between these two methods, which is more notice- 
able for control parameters (T, e, 9 and to) of order of 1. This is expected since 
the analytical approximation to the eigenvalue problem should be accurate for con- 
trol parameters significantly less than 1 (see Appendix 16.11 ). In Figs. [T]l3]we also 
plot the evolutionary timescales obtained from the order-of-magnitude estimates in 
Sect. 14. II (open circles) [see Eqs. (l60l ). (|6TI ) and (l62l )1 with numerical coefficients 
inserted by hand to ensure reasonable agreement with the other, more precise de- 
terminations, 

h = lA/ilnf, (68) 



h = 0Jl[il7:fe + e 



1 



1 



1 



7.5rW^)^e 2.59 



+ 0.Sl{l7rrZ] 



(69) 



h = 2.3 



-I- 



[(Z7r)2? + e] 
1 / 1 



-1-1 



1 



7.5Y Xil^r^ 2.59 



+ 0.8I(Z;r)^w 



(70) 



The fact that this agreement can be reached using coefficients that do not differ 
significantly from unity (see Eqs. ( [68] ). ( |69l ) and dTOl )) corroborates the order-of- 
magnitude estimates in the framework of the linear-mode regime. 
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Figure 1 : Characteristic time 7i on which the system achieves magnetohydrostatic 
quasi-equilibrium, as a function of the dimensionless control parameters of the 
different physical processes (i.e., T [magnetic field strength], e [collision rate], 
9 [weak interaction rate] and to [ohmic resistivity]). The time variable t is normal- 
ized as 7 = t/ts where = acf l{nQ„knB) and a is an artificial parameter controlling 
how rapidly the system reaches magnetohydrostatic quasi-equilibrium and knB is 
a property of the background system. We show a comparison of the characteristic 
time obtained from the order-of-magnitude estimate of Eq. (1681 ) (open circles) with 
those obtained from the numerical solution of the eigenvalue problem of Eq. (1 100 1 ) 
(black diamonds) and from the analytical approximation of Eq. (11231 ) (open trian- 
gles), (a) 7i as a function of Y for ? ^ 0.02, ^ = 0.05, and w = 0.0001. (b) h as 
a function of ?for Y = 0.01, 'B = 0.05, and oj = 0.0001. (c) h as a function of 9 
for ? = 0.02, Y = 0.01, and to = 0.0001. (d) h as a function of w for e = 0.02, 
Y = 0.01, and ^ = 0.05. 
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Figure 2: Same as Fig. [T]but for the characteristic time ?2 on which the particle 
densities reach their diffusive quasi-equilibrium. We show a comparison of the 
characteristic time obtained from the order-of-magnitude estimate Eq. (l69l ) (open 
circles) with that obtained from the numerical solution of the eigenvalue problem 
of Eq. (1 100 1 ) (black diamonds) and from the analytical approximation of Eq. (11241 ) 
(open triangles). 
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Figure 3: Same as Fig. [T] but for the characteristic time on which the magnetic 
field evolves. We show a comparison of the characteristic time obtained from the 
order-of-magnitude estimate Eq. (T/Ol) (open circles) with that obtained from the 
numerical solution of the eigenvalue problem of Eq. (1 1001) (black diamonds) and 
from the analytical approximation of Eq. (|124b (open triangles). 
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4.3 Finite difference evolution of a linear magnetic-field perturbation 

We developed a numerical code based on a finite-difference algorithm that solves 
the nonlinear set of Eqs. 02] ). ( [33] ). and ( [34l ) (see Appendix I6.2I ). In this section, 
we use that code to simulate the evolution of a linear magnetic-field perturbation 
of the form given by Eq. ( [64l) . as well as the associated particle-density evolution. 
In this task, we aim to achieve two main goals: 

1 . To verify the quasi-equilibrium states that occur during the evolution. 

2. To compare the linear normal-mode solution given by Eqs. (|65]) . (l66l ). and (l67] ) 
with the solution obtained by solving numerically the set of Eqs. ([32]), ([33]) and 
([34])(see Appendix 16.21 ), to measure the precision of the numerical 
code. 

In our numerical simulation, we set the initial condition to correspond to a mag- 
netic field of the form given by Eq. ([64]) . with 6B.(x,0) = Acos{lnx), A = -0.010 
and 1-2. We set the initial particle density perturbations to be 5nc(x,0) = 
Sn„(x,0) = 0. In what follows, we neglect the ohmic diffusion, therefore, we 
expect that t\ = Imhs^ ~h - 'tsuc^ Is = 1b- We use the parameters: Hj - 0, 
T = 0.20, 6 - 0.\0,'e - 0.010. We note we are in a regime in which ambipolar 
diffusion occurs more rapidly than the weak interaction processes Qdrag ~ 2.5 and 
'tweak ~ 10). From a numerical solution of the eigen- value problem [see Eq. (11001) 1. 
we obtain 1\ = 0.028, t2 = 1.2, and ?3 = 39.0. In Fig. [4] we show the system's 
evolution during the first timescale li , which is the timescale over which the system 
can reaches magnetohydrostatic quasi-equilibrium. We label the different instants 
of this evolution with progressive numbers as described in the figure. This mag- 
netic perturbation perturbs the hydrostatic quasi-equilibrium of the homogeneous 
background star; thus particles are therefore compelled to move in order to com- 
pensate for this imbalance, and reach magnetohydrostatic quasi-equilibrium in a 
short timescale of the order of 7i . In Fig. [4] (a) we observe the evolution of the 
magnetic-field perturbation during this short timescale, and that it has not evolved 
significantly. This is expected since the magnetic pressure is small compared with 
the background-degeneracy pressure, and can therefore induce only a small rela- 
tive density perturbation in the particles (see Figs. [4] (b) and [4] (c)). We can also 
see in Figs. [4] (b) and[4](c) that during this short timescale dndriQc ~ 3n„/no„, 
as expected from Eq. ([38] ). In Fig. [5] we show the system's evolution during the 
second timescale I2, on which the diffusive quasi-equilibrium of the charged parti- 
cles is attained. At the end of this timescale, 5?i„/?ion ~ 0, and dndriQc has grown 
significantly to balance the magnetic pressure. Towards the end of this timescale, 
the separation between the Unes is narrower, showing that the charged particles 
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have reached their diffusive quasi-equilibrium. In Fig. [6l we show the system's 
evolution during the third timescale 73, on which the chemical quasi-equilibrium 
is restored and the magnetic perturbation decays. At the end of this timescale, 
the density perturbations as well as the magnetic-field perturbation have decreased 
substantially. 

On the other hand, we note that a linear magnetic-field perturbation can be 
written in general, as a superposition of the linear normal modes given by Eq. (|67] ) 
as 

3 

5B^(x, ~t) = Y^ ^massjjim) cosilnx) exp(-?7,„?). (71) 

m=l 

where rj^ = l/7„, and the coefficients D,„ are calculated from the initial condi- 
tions 6B^(x, 0), Sndx, 0), and dn„(x, 0). For the initial conditions of this simulation 
and from the numerical solution of Eq. (llOOl ). we derive that a^^ (jji) = -0.0076, 

(^2) = 0-10. «5B^('73) ^ 0-22, Di = 0.023, D2 - -0.016 and D3 - -0.037. We 
Note that, for 7 » 72, 

5B^(x, 7) K D^a^^ (773) cos(/7r7) exp(-7737). (72) 

On the other hand, the total magnetic energy associated with this perturbation 
in physical units is 

6EB{t)^ f ^^§^dx. (73) 
.in 6n 



Thus, in the asymptotic limit of Eq. 

log{6EB{t)) « log{6E'^'m - 2 log(^')-. (74) 



t3 



where SEfiO) = (DsaseXm)) d/(l67T) is the initial energy contained in the third 
mode. To measure the precision of our numerical code, in Fig. |7] we compare the 
evolution of the magnetic-energy perturbation obtained from the simulation with 
the asymptotic limit of Eq. (1741 ). In Fig. |7j we see that there is good agreement 
between our numerical result and the asymptotic limit of Eq. d74l ). 



4.4 Finite difference evolution of a non-linear magnetic field 

Using our numerical code, we simulate the evolution of a non-linear magnetic field, 
by applying the set of eqs. (I32l ). (1331 ). and (l34l ) without any linearization of B^. The 
main aims of this section are: 
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Figure 4: Evolution of a magnetic field of the form B,{x, t) = Be + 6B^{x, t), where 
Be is a constant and 6B^ is a small perturbation (i.e., I^B^I «; Be). The associated 
charged-particle and neutrons density perturbations around a homogeneous back- 
ground star are dndx, t) and 5nn{x, t) respectively, with the neutron and charged- 
particle number densities in the background defined by «()„, hqc, and d the total 
length of the system. This figure shows the evolution during the first timescale ti 
on which the system reaches the magnetohydrostatic quasi-equilibrium. The time 
variable is normalized as in Figs. [B [2] and [3] We have used in this simulation 
riQn - 9.8 X 10^'cm"^ and noc/non = 4.0 x 10^^. We set the control parameters 
of the different physical processes to: 0, T = 0.2, 9 - 0.1, and e - 0.01, 
which gives 7i = 0.028. (a) Magnetic field perturbation at different instants labeled 
with progressive numbers starting with (1) for the initial condition at 7 = 0. The 
initial magnetic perturbation is dB^{x,0) = BcAcos(lnx/d), with A = -0.010 and 
1 = 2. while the initial particle density perturbations are 6nc{x,0) = 0) = 0. 

The other instants are(2): 7 = 0.014, (3): 7 = 0.028. (b) Charged-particle density 
perturbations at the same instants as in panel (a), (c) Neutron density perturbations 
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Figure 5: Same parameters and normalization conventions than Fig.|4]but showing 
the evolution during the second timescale ?2 on which the particles densities evolve. 
For these parameters t2 = 1-2. (a) Magnetic field perturbation at different instants. 
The instant labels are: (3): t = 0.028, (4): 7 = 2.4, (5): 7 = 3.8, (6): 7 = 4.2, (7): 
7 = 4.7. (b) Charged-particle density perturbations at the same instants as in panel 
(a), (c) Neutron density perturbations at the same instants as in panel (a). 
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Figure 6: Same parameters and normalization conventions than Fig.|4]but showing 
the evolution during the third timescale 1^ on which the magnetic evolves signifi- 
cantly. For these parameters ?3 = 39.0. (a) Magnetic field perturbation at different 
instants. The instant labels are: (7): t = 4.7, (8): t ^ 27.0, (9): 7 = 116.0, 
(b) Charged-particle density perturbations at the same instants as in panel (a), (c) 
Neutron density perturbations at the same instants as in panel (a). 
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Figure 7: Decay of the total magnetic energy contained in the perturbation of Fig.|4] 
(a), Fig.|5](a), and Fig.[6](a) SEsit) = (l/S/r) j^{6B^fdx (points). We compare this 
solution with the normal mode solution in the asymptotic limit of Eq. (1741 ) (full 
line). The quantity 6E'-g \0) is the initial energy contained in the third mode whose 
decay time is ti,. The time variable is normalized respect to the decay time of the 
second mode tj. For the parameters of this simulation t3/t2 = 32.5. 
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1. To compare the estimate obtained using Eq. ( [59l ) with the resuh of the nu- 
merical code in the non-linear regime. 

2. To verify the conservation laws given by Eq. (I20l) for the magnetic flux and 
baryon number, which provide a measure of the right performance of our 
numerical code. 

Defining the initial conditions of our simulation, we set the particle densities, 
Sndx, 0) = SEnix, 0) = 0, for an initial Gaussian magnetic field profile given by 



We normalize the magnetic field in Eq. (1751 ) to its maximum value B"'"-''. The 
parameter s in Eq. dTSl) is related to its standard deviation p by s-^ = d^/2p^, 
and xq = xo/d is the center of the Gaussian function. The characteristic spatial 
variation in the Gaussian function can be estimated by the width w of the Gaus- 
sian function at its half-height, which is given by w w 2.355/?; we can therefore 
approximate the characteristic length over which the magnetic force varies to be 
L ~ 2355d/^/2s. In our analysis, we again neglect ohmic diffusion. To 
estimate the characteristic decay time of the magnetic field in Eq. (T/Sl ). we use as 
a first approximation the characteristic length of the magnetic force in Eq. ( [59l ) as 
the width of the initial gaussian, therefore. 



In our numerical analysis, we assume that s - 20.0, xo - 0.5, To - 0,e = 0.01, 
9 = 0.1, and Y = 0.2. If we substitute these parameters in Eq. (1761 ). we obtain a 
characteristic decay time Ib « 53.5. In Fig.[8](a), we see the evolution of the initial 
magnetic field given by Eq. (1751) . and that at times close to Ib the maximum of 
the magnetic field has decayed to 60% of its initial amplitude. On the other hand. 
Fig. [8] (b) shows the evolution of the charged-particle density perturbations, and 
Fig- [11(c) shows the evolution of the neutron density perturbations. 

We verified the conservation law for the dimensionless baryon number NbO) = 
Nob + SNbO), where Nqb - 1-04 was the backgound baryon number. According to 
the initial conditions SNBiO) = 0, which implies that the conservation law requires 
NbO) = A^s(O) = Nob- The calculated values of the baryon number in the simula- 
tion time steps show a maximum percentage error of 2.0 x 10"^ % with respect to 
Nob- Regarding the magnetic flux, we found exact conservation with respect to the 
initial magnetic flux 0^(0) down to 7-digit precision. 



0) ^ exp{-s^{x - xof). 



(75) 




(76) 
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Figure 8: Evolution of an initial Gaussian magnetic field profile given by B^{x, 0) = 
^max Qxp{-{20/d)^{x-{d/2))^) as well the associated particle density perturbations. 
We have used the same background parameters and time normalization convention 
as in Fig. |4l The initial condition in the particle densities is such that they are 
uniformly distributed in the background, thus dn„{x,0) = dnc{x,0) = 0. We set 
the parameters controling the different physical processes as: co = 0, e = 0.01, 
9 = OA, and T = 0.2. The characteristic decay time over which the magnetic field 
decays significantly is Ib ~ 53.5. (a) Evolution of the magnetic field for different 
instants labeled with progresive numbers, (1) : 7 = 0. (2) : 7 = 54, (3) : 7 = 180, 
(4) : 7 ^ 500, (5) : 7 ^ 1000. (b) feVolution of the charged-particle density 
perturbations for the same instants as in panel (a), (c) Evolution of the neutron 
density perturbations for the same instants as in panel (a). 





5 Summary and Conclusions 



We have studied the long-term evolution of the magnetic field and the densities of 
neutrons and charged particles (protons and electrons in the interior of a neutron 
star, using a multi-fluid model with a simphfied geometry in which the magnetic 
field points in one Cartesian direction and varies along an orthogonal direction. 
We found a set of three non-hnear partial difi'erential equations of first order in 
time describing this evolution, and we analyzed them in three different ways: (i) 
estimating evolutionary timescales directly from the equations, guided by physical 
intuition; (ii) a normal-mode analysis of the equations in the limit of a nearly uni- 
form system; and (iii) a finite-difference numerical integration of the equations. In 
this section, we summarize the main results of each one of these approaches and 
present the main conclusions of this work. 

5.1 Evolutionary timescales 

The three partial differential equations of our model constitute a dynamical system 
of three independent variables (magnetic-field, charged-particle, and neutron den- 
sity perturbations). We identified three characteristic evolutionary timescales on 
which the system approaches successive quasi-equiUbrium states. In the follow- 
ing, we summarize the basic physical processes governing the evolution and the 
estimates of each timescale. 

5.1.1 Timescale to achieve magnetohydrostatic quasi-equilibrium 

During the early stages of a neutron star's life, the Lorentz force moves the bulk 
stellar fluid, inducing perturbations in its pressure. In this short time, all parti- 
cles move together as a single fluid with the same bulk velocity. There is also 
not enough time for weak interactions between particles (beta decays) to operate, 
so the composition is frozen. The system evolves until it reaches a magnetohy- 
drodystatic quasi-equilibrium state, in which the Lorentz and the fluid forces are 
close to balancing. Neutron stars with no rotation or convection reach this quasi- 
equilibrium state in a short time, not much longer than the Alfven time, which for 
typical magnetar core parameters scales as 

tAlfven = 5.7 X lO'^ R(, B\l s, (77) 

where = R/{10^ cm) denotes the radius of the star in units of 10^ cm and 

Bi5 = B^/{\0^^ G) the magnetic field in units of lO'^ G. On timescales far longer 
than the Alfven time, relative movements between the different species of particles 
in the star as well as weak interactions between them become relevant (GR-92). 
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Therefore, a description of the system during these stages requires a multi-fluid 
theory. In this paper, we develop a model to study the decay of the magnetic field 
induced by these long-term mechanisms. The Alfven time is far shorter than the 
timescales of these long-term processes (see Sects. 15.1.21 and 15.1.31 ). which are 
of the order of years or much more. A numerical code simulating the evolution 
during the Alfven timescale would require a time step many orders of magnitude 
shorter than that required to simulate the long-term evolution in a computational 
time not prohibitively long. In our model, we overcome this difficulty by replacing 
the short-term dynamics by a ficticious friction force term acting on the neutrons 
(the most abundant species). Its strength is controlled by an artificial parameter a, 
chosen so that the timescale to reach magnetohydrostatic quasi-equilibrium is long 
enough for the numerical code to be able to deal with it (and therefore much longer 
than the Alfven time), but shorter than the timescales of the long-term processes 
that we are interested in modeling. 

5.1.2 Timescale for the evolution of the particle densities 

On timescales far longer than required to reach magnetohydrostatic quasi-equilibrium, 
the neutrons and charged particles move relative to each other with different veloc- 
ities and are affected by coUisional drag. Also, weak interactions (beta decays) 
convert particles from one species into another, erasing departures from chemical 
quasi-equilibrium among neutrons, protons, and electrons caused by the induced 
particle density perturbations. Both processes contribute to the evolution of the 
particle densities. In estimating the timescale for this evolution, it is useful to 
distinguish the regimes in which each one of these processes is dominant. If the 
collisions between the charged particles and the neutrons are rare and the weak in- 
teractions are slow, the charged particles move easily relative to the neutrons (am- 
bipolar diffusion) allowing the particles to reach a diffusive quasi-equilibrium state 
(in which gravitational, electromagnetic, and pressure forces are closely balanced) 
but staying avoid chemical quasi-equilibrium. The diffusive quasi-equilibrium 
is reached in a timescale controlled by the collision rate between neutrons and 
charged particles given by 

tdrag- "^.5 X 10-' L^T^yr, (78) 

where Lj = L/{\0^ cm) and = 7/(10^ K). 

In the opposite regime, when the collisions between charged particles and neu- 
trons are frequent, and the weak interaction rate is high, the relative diffusion of 
particles will be impeded by the coUisional drag, but the system reaches the chem- 
ical quasi-equilibrium in a timescale controlled by weak interactions of magnitude 
(assuming modified Urea reactions) 
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W ~ 4.3 X 10^ yr. 



(79) 



5.1.3 Timescale for magnetic field evolution 

The Ohmic dissipation promotes the decay of the magnetic field in a timescale 

tohmic ~ 1.4 X 10^' LjT-^yr. (80) 

This extremely long time implies that Ohmic dissipation is not very effective in 
producing magnetic-field decay in neutron stars. A more rapid decay of the mag- 
netic field can occur through ambipolar diffusion, i.e., a joint drift of the magnetic 
field and the charged particles relative to the neutrons, in which magnetic energy 
becomes dissipated by collisions. This diffusion process induces local density per- 
turbations that deviate from chemical quasi-equilibrium. Weak interactions tend to 
erase the chemical imbalance by converting particles of one species into the other. 
During these conversions, neutrinos, and antineutrinos remove part of the mag- 
netic energy. In estimating the timescale of magnetic field decay induced by these 
processes, it is useful to distinguish between the two regimes we discussed above. 

If ambipolar diffusion occurs far more rapidly than the weak interaction pro- 
cess, i.e., tiirag ^ t„eak, diffusivc quasi-cquilibrium is quickly reached, but chem- 
ical quasi-equilibrium is not. Thus, the rate at which chemical quasi-equilibrium 
is restored by weak interactions determines the characteristic timescale over which 
the magnetic field evolves, 

tB ~ 1.7 X 10"^ B-^ T^^ yr. (81) 

This timescale is the corresponding one to reach chemical quasi-equilibrium t^eak 
but amplified by a factor of the order of the ratio of the charged fluid pressure to 
the magnetic pressure. This is expected since the Lorentz force drives the am- 
bipolar diffusion that prevents weak interactions from restoring chemical quasi- 
equilibrium. 

In the opposite regime, collisions are frequent and weak interactions are fast, 
i.e. tweak ^ tdrag, SO the systcm reaches chemical quasi-equilibrium long before 
diffusive quasi-equilibrium. Therefore, the evolution of the magnetic field is lim- 
ited by the collisions that control the ambipolar diffusion, on a timescale of 

tB ~ 1.8 X IQ"' B^jLlTlyr. (82) 

Since the magnetic field sustains the diffusive imbalance and drives the diffusion, 
the timescale given in Eq. (I82l ) is similar to that taken by particles to achieve dif- 
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fusive quasi-equilibrium but, as above, amplified by a factor that depends on the 
magnetic field strength. 

In the limit of non-interacting protons and neutrons studied by GR-92, Eqs. (ISTI ) 
and (|82l ) correspond to the estimate of the timescale for magnetic field evolution 
due to the ambipolar diffusion process of Eq. (35) in GR-92. We note that GR- 
92 considered the neutrons to be a static species in diffusive quasi-equilibrium. 
From this correspondence, we confirm that the motion of the neutrons, which was 
included in our model, plays no important role in the one-dimensional evolution. 

5.2 Normal-mode solutions 

In the spirit of finding an analytical solution to the set of non-linear partial differ- 
ential equations that describes the evolution of the system, we completed a linear 
perturbation analysis of the equations in the limit of a nearly uniform system and 
tried a normal mode solution (Sect. 14.21 ). We found three exponentially-decaying 
modes that could be interpreted as the episodes of successive relaxation to the 
quasi-equilibrium states identified in the previous physical analysis (see Sects. 14.11 
and 15.11 ). We found analytical and numerical solutions for the decay times of these 
modes (see Appendix l6.1l ) that agree with the estimated formulae for the evolution- 
ary timescales of the general non-linear system (see Sects. 14. 1[ 15.11 and Figs. I1I2[ 
and El). 

5.3 Numerical integration of the partial differential equations 

We constructed a finite-difference numerical code to solve the full system of non- 
linear partial differential equations. This was tested against the normal-mode solu- 
tions (Fig. IT]) and by verifying the known conservation laws of magnetic flux and 
baryon number in highly non-linear situations (Fig. [8]). Applying this numerical 
code, we propose to consider, in the future, conditions in which the background is 
non-homogeneous, as well as non-linear situations. 

5.4 Conclusions 

We have established a general multifluid formalism in which it is possible to study 
the magnetic field decay processes in neutron stars identified by GR-92, which 
are likely to represent the basis of the magnetar phenomenon. As a first step, we 
have focused on a simplified geometry in which the magnetic field points in one 
Cartesian direction and varies in an orthogonal direction to it. We estimated the 
timescales of the relevant processes, and using numerical simulations we followed 
the temporal evolution of some simple magnetic-field configurations. The present 
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work is far from exhausting the possibilities of this formalism, which should be 
applied to more realistic situations, including neutron stars with spherical symme- 
try, non-uniform background density and three-dimensional magnetic fields with 
components that depend on two or three spatial coordinates. Further steps might 
consider convective motions and additional species of particles as well as the ef- 
fects of superfluidity and superconductivity. 



6 Appendix 



6.1 Solution to the linear system 

The set of Eqs. (I32l) . (I33l ). and (l34l ) is, in general, non-linear with respect to the 
magnetic field variable B^. To find a linear solution to this set, we linearized the 
equations respect to the magnetic field, as described in Sec. 14.21 Assuming that 
the properties of the background star are homogeneous, that is quantities with 
sub-index zero do not depend on position, we obtain the following linear set of 
equations: 



^ = Do-^ + Wo Q, 



where we define the solution vector to be 



with the coupling matrices 



Q 



SB, > 



(83) 



(84) 



where 



Dn 



d3i 



dn di3 " 
dii d23 
d32 d33 , 



dn = (1 ?i)[l yS(l + 



\ + n 
e 



in 



k I + n 



)], 
)], 



di3 = 2{\+n)[l + 



I + n 



], 



«[l+y6(l + -)], 

n 



(85) 

(86) 
(87) 
(88) 
(89) 



36 



d22 = kn[l + f (1 + -)], 
k n 



d23 = 2n[l + -], 
n 



dji = Td2l, dT,2 = Td22, <5?33 = Ci> + Td23 



(90) 
(91) 
(92) 



and 



where 



Wn 



Wll Wl2 Wl3 
W21 1^22 VV23 
^^ ^31 H'32 ^33 J 



Wn - W\2 - Wi3 = H'23 = W31 - H'32 - Wt,3 - 0, 



(93) 



(94) 



W21 ^ -n9{- - 1), 



(95) 



^22 = -n6{ k). 

n 



(96) 



To simplify the notation, we define cJ = a),e = e, T = T, 9 = B,rQ = r, hqc = n, 
knc = k, g = nkcc, (3 = nkps- We find normal-mode solutions to this linear system 
of the form 



= exp(-77„,7)cos(/7r;c)Am, 



(97) 



where m = 1,2,3 labels each mode, and we order the decay times in increasing 
order. The decay time for each mode is 



and the amplitudes are 



A - 



(98) 



(99) 



where is normalized with respect to while agngc both normalized with 



respect to 5ns - B^/(SnnonknB)- If we substitute Eq. (|97| ) into the linear matrix 
system Eq. (|83] ). we obtain the eigenvalue problem 
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Co Am — 77^A,„, 



that needs to be solved, where 



Co - {InYBo - Wo. 



(100) 



(101) 



6.2 Analytical solution for the decay times 

To obtain analytical expressions for the decay times of the normal modes [see 
Eq. (|98])1. we write the matrix of Eq. (llOll l as 



Co = 



Cll Cl2 Ci3 

UCu +X21 UCn +X22 UCu +X23 

{ fcn +^31 fcn +A'32 fcu +^^33 } 



where we define n* = \ + n, u = n/n*, e* = e/n*, f = uT, and 

X2i=ilnf/ie* +ne(^-iy 

^2 * '^/?^23 



X32 ^ Z , X33 = TY23 + (ot) 0). 



(102) 

(103) 
(104) 
(105) 
(106) 
(107) 
(108) 
(109) 



We also note that all;^';; «: 1 whereas Cjj ~ 1. In the limit = 0, the eigenvalues 
are ?7i = cn + ucn + fc\3, 772 = ??3 = 0. In the general case, we can write 



771 = cn+ucn+ fc\3 + 9 



= {lnfn[\+fi + n{2r + k) 
+e*(J3 + n{2r+p))]+(p, 



(1 10) 
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and treat 7/2, t]3, (fi ~ Oixij) as small quantities, which we derive below. We calcu- 
late the characteristic polynomial of the matrix Eq. ( 11021 ) and compare terms of the 
same order in this polynomial with those in the equation 

(m - - mm - Ji) - 0. (iii) 

After neglecting quadratic and cubic terms in the small quantities 772, ?73, and(p, we 
find 

V2 + m = F, (112) 

mm = G, (113) 

'P = (X22+X33)-F, (114) 



where 



cnF\ - C12F2 - cuF3 



il7T)^n*[l +P + niir + k) + e*{p + ?i(2Y + p))] ' 

with the definitions, 



and 



with. 



(115) 



F\ =X2i +X2-i, (116) 

F2 =X2l + /Y23 - "A'33, (117) 
^^3 = UX32 +X'i\ - fX22, (118) 



Q_ ciiGi +C12G2 + C13G3 

" (/;r)V[l +yS - «(2Y + /:) + e*06 - «(2T +p))] ' 



G 1 ^ X22X?,?, - X2-iXz2 , ( 1 20) 

Gi ^X23X3i -X2\Xi>-i^ (121) 

G3 ^^21^32 -X22XT,\- (122) 
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Thus, the eigenvalues are given by 

7/1 - (InfnU +15 + «(2Y + k) + e{{i + «(2T +p))] + if. 



(123) 



F± iF^- AG 

?72, m = 2 ■ ^^^^^ 

Equations. (11231 ) and (11241 ) give approximate analytical expressions for the decay 
times of the normal modes [see Eq. ( |98l) 1 as functions of the dimensionless pa- 
rameters, T„,(e, T, 9, a>). Equations. (11231 ) and (11241 ) are useful because from them 
we can understand in an analytical way the relative importance of the different 
physical processes for the behavior of the decay times. The linear decay times fm 
correspond to the estimated characteristic timescales t,„ in the general non-linear 
regime (see Eqs. (l60l ). (|6TI ). and (l62l)). 



6.3 Numerical Method 

To solve the dimensionless system of Eqs. (l32l) . (l33l) . and (l34l ) numerically we dis- 
cretize them using an FTCS (Forward Time Centered Space) first order scheme. We 
define the spatial and temporal grid by xj = jAx, ~t„ = nAt where 7 = 0, 1, 2, j,„ax 
and n = 0, 1,2, ...,n,„ax- In this method, we evaluate the spatial derivatives of a 
quantity A(x,7) at a half point between consecutive grid points, namely, 

dA\ A._,i,^ ^^^^^ 



(126) 



dx I J Ax 
dA\ 



- ' (127) 



dx lj_i/2 

Then we interpolate it to the grid using the known grid quantities, 

A" , + A" 

Ah/2 = ^V^' ^^^^^ 

A" +A" 
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